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INVARIANT DOMAINS AND FIRST-ORDER 
CONTINUOUS FINITE ELEMENT APPROXIMATION 
FOR HYPERBOLIC SYSTEMS* 

JEAN-LUC GUERMONDt AND BOJAN POPOVt 

Abstract. We propose a numerical method to solve general hyperbolic systems in any space 
dimension using forward Euler time stepping and continuous finite elements on non-uniform grids. 
The properties of the method are based on the introduction of an artificial dissipation that is defined 
so that any convex invariant sets containing the initial data is an invariant domain for the method. 
The invariant domain property is proved for any hyperbolic system provided a CFL condition holds. 
The solution is also shown to satisfy a discrete entropy inequality for every admissible entropy of the 
system. The method is formally first-order accurate in space and can be made high-order in time by 
using Strong Stability Preserving algorithms. This technique extends to continuous finite elements 
the work of [Hoff(1979)[|Hoff(1985)] , and [Frid(2001)] . 
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1. Introduction. The objective of this paper is to investigate a first-order ap¬ 
proximation technique for nonlinear hyperbolic systems using continuous finite el¬ 
ements and explicit time stepping on non-uniform meshes. Consider the following 
hyperbolic system in conservation form 


dtu+ \/-f{u) = 0, for (a;, t) e 
u(x,0) = Uq(x), fora^eM'^. 


where the dependent variable u takes values in K"* and the flux f takes values in 
In this paper u is considered as a column vector u = (iti,..., UmY■ The flux 
is a matrix with entries (tt), l<j<d and V-/ is a column vector with 

entries (V-/)i = ^or any n = (ni... ,nd)^ G we denote f{u)-n 

the column vector with entries where i € The unit sphere 

in centered at 0 is denoted by 1). 

To simplify questions regarding boundary conditions, we assume that either pe¬ 
riodic boundary conditions are enforced, or the initial data is compactly supported 
or constant outside a compact set. In both cases we denote by D the spatial domain 
where the approximation is constructed. The domain D is the d-torus in the case of 
periodic boundary conditions. In the case of the Cauchy problem, D is a compact, 
polygonal portion of large enough so that the domain of influence of Uq is always 
included in D over the entire duration of the simulation. 

The method that we propose is explicit in time and uses continuous finite elements 
on non-uniform grids in any space dimension. The algorithm is described in §3.2[ see 
(3.5) with definitions (3.4)-(3.8l-(3.13|. It is a somewhat loose adaptation of the 
non-staggered Lax-Friedrichs scheme to continuous finite elements. The key results 
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of the paper are Theorem |4.1| and Theorem |4.2| It is shown in Theorem |4.1| that the 
proposed scheme preserves all the convex invariant sets as dehned in Definition |2.3 
and it is shown in Theorem |4.2| that the approximate solution satisfies a discrete 
entropy inequality for every entropy pair of the hyperbolic system. Similar results 
have been established for various hnite volumes schemes by |Hoff(197^ Hoff(198^ , 
Perthame and Shu(1996)], |Frid(200'T)] for the compressible Euler equations and the 


p-system. Our scheme has no restriction on the nature of the hyperbolic system, 
besides the speed of propagation being finite. To the best of our knowledge, we are 
not aware of any similar scheme in the continuous finite element literature. 

The paper is organized as follows. The notions of invariant sets and invariant 
domains with various examples and other preliminaries are introduced in Section 
The method is introduces in Section Stability properties of the algorithm are 
analyzed in Section Numerical illustrations and comparisons with existing first- 
order methods are presented in Section 

2. Preliminaries. The objective of this section is to introduce notation and 
preliminary results that will be useful in the rest of the paper. We mostly use 
the notation and the terminology of [Chueh et al.(1977)Chueh, Conley, and Smoller 


Hoff(197^ Hoff(1985)| Frid(200T)] . The reader who is familiar with the notions of 


invariant domains and Riemann problems may skip this section and go directly to 
although the reader should be aware that our definitions of invariant sets and domains 


are slightly different from those of Chueh et al.(1977)Chueh, Conley, and Smoller 
Hoff(1979)| [Hoff(1985)t[Frid(2001)| . 


2.1. Riemann problem. We assume that (1.11 is such that there is a clear 
notion for the solution of the Riemann problem. That is to say there exists an 
(nonempty) admissible set A C K"* such that for any pair of states [ul^ur) £ Ay, A 
and any unit vector n G 1), the following one-dimensional Riemann problem 


(2.1) dfU + dx(f(u)-n) = 0, (a;,<) € MxIR+, tt(a;,0) 


Ur, if X < 0 
Ur, if X > 0, 


has a unique (physical) solution, which we henceforth denote u(n,UR,UR). 

The theory of the Riemann problem for general nonlinear hyperbolic systems with 
data far apart is an open problem. Moreover, it is unrealistic to expect a general the¬ 
ory for any system with arbitrary initial data. However, when the system is strictly 
hyperbolic with smooth flux and all the characteristic fields are either genuinely non¬ 
linear or linearly degenerate, it is possible to show that there exists d > 0 such that 
the Riemann problem has a unique self-similar weak solution in Lax’s form for any 
initial data such that \\uR — UR\\t 2 < 5, see [Lax(1957)] and |Bressan(2000) Thm 5.3]. 
In particular there are 2m numbers 


( 2 . 2 ) 


Ar < A+ < < A+ < ... < A- < A+ 


dehning up to 2m -I- 1 sectors (some could be empty) in the (x, t) plane: 


(2-3) -g(-oo,A;^), - e (A;^ , A]*"),..., -€(Aj,^,A+), - e (A+,oo). 


The Riemann solution is ur in the sector | G (—oo, Aj^ ) and ur in the last sector 
f e (A+ ,oo). The solution in the other sectors is either a constant state or an 
expansion, see [Bressan(2000)| Chap. 5]. The sector A^t < x < A+t, 0 < t, is 
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henceforth referred to as the Riemann fan. The key result that we are going to use 
is that there is a maximum speed of propagation A„iax(’^; 14^;) := max(|A^|, |A+|) 

such that for t > 0 we have 


(2.4) 


u{x t) - if a; < -iAmax(«,ML,Mii:) 

[Mij, if a; > tAmax(«-,ML,Mii)- 


Actually, even if the above structure of the Riemann solution is not available or 
valid, we henceforth make the following assumption: 

The unique solution of (13 has a finite speed of propagation for any n, 
i.e., there is Aniax(^i, “l, Mfl) such that (2.4) holds. 


(2.5) 


For instance, this is the case for strictly hyperbolic systems that may have char¬ 
acteristic families that are either not genuinely nonlinear or not linearly degener¬ 
ate, see e.g., |Liu(1975')| Thm.1.2] and [Dafermos(2000)[ Thm. 9.5.1]. We refer to 
|Osher(1983)[ Thm. 1] for the theory of the Riemann problem for scalar conservation 
equations with nonconvex fluxes. In the case of general hyperbolic systems, we re¬ 
fer to [Bianchini and Bressan(2005) Section 14] for characterizations of the Riemann 
solution using viscosity regularization. We also refer to [Young(2002)[ Thm. 2] for 
the theory of the Riemann problem for the p-system with arbitrary data (i.e., with 
possible formation of vacuum). 


The following elementary result is an important, well-known, consequence of (2.4) 


i.e., the Riemann solution is equal to ul for x € (—oo, Aj^ t) and equal ur for x S 
(A+t,oo): 

Lemma 2.1. Let UL,Ufi G A, let u{n,UL,Uf{) be the Riemann solution to 

let u(t, n, Ul, ub) '■= Pi u{n, ul,ub)(x, t) dx and assume that t Xmuxin, ul,ub) < 
2 

then 

( 2 . 6 ) u{t,n,UL,UB) = ]^[ul + Up.) - t{f{uB)-n - f{uL)-n). 


If the system (l.I) has an entropy pair (rj,q), and if the Riemann solution is 


defined to be entropy satisfying, i.e., if the following holds 

(2.7) dtf]{u{n,UL,UB)) + d,,{q{u{n,UL,UB))-n) <0, 

in some appropriate sense (distribution sense, measure sense, etc.), then we have the 
following additional result. 

Lemma 2.2. Let {r],q) be an entropy pair for 0) and assume that \2.1\ holds. 


Let Up, Up G A and let u{n,up,upi) be the Riemann solution to (2.1). Assume that 
^ -^max {n,up,upi) < A, Then 

( 2 . 8 ) f]{u{t,n,UL,UB)) < \{r]{uL) + q{un))-t{q{un)-n-q{up)-n). 


Proof. Under the CFL assumption t\m_a,x[n,up,upi) < h the inequality ( |2.7| ) 
implies that 


(2.9) 


■q{u{n,UL,UR)){x,t) dx < \{q{uL) + r]{un)) - t{q{un)-n - q{uL)-n). 


Jensen’s inequality p(M(t, n, ml, it^j)) < pi q{u{n,up,upi)(x,t)) dx then implies the 

2 

desired result. □ 
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2.2. Invariant sets and domains. We introduce in this section the notions 
of invariant sets and invariant domains. Our definitions are slightly different from 
those in [Chueh et al.(1977)Chueh, Conley, and Smoller[ Hoff(1985), Smoller(1983) 


Frid(2001)| . We will associate invariant sets only with solutions of Riemann problems 
and dehne invariant domains only for an approximation process. 

Definition 2.3 (Invariant set). We say that a set A C A C K™ is invariant for 
ifforanypair{ui^,Uii)£AxA, any unit vector n G ^(0,1), andanyt>0, 


the average of the entropy solution of the Riemann problem (2.1) over the Riemann 
fan, say, , ^ M(n, ml, t) dx, remains in A. 

t ^ y\ A 2^ ) A 2 t 

Note that, the above definition implies that given t > 0 and any interval I such 
that {Xft, A+t) C I, we have that j u{n,UL, ub){x, t) dx G A. Note also that most 
of the time expansion waves and shocks are not invariant sets. 

We now introduce the notion of invariant domain for an approximation process. 
Let X/j C be a finite-dimensional approximation space and let Sh '■ Xh 9 

Uh I— Sh{uh) G Xfi be a discrete process over X^. Henceforth we abuse the language 
by saying that a member of X^, say Uh, is in the set A C K"* when actually we mean 
that {uh{x) I a; G M} C H. 

Definition 2.4 (Invariant domain). A convex invariant set A C A C M"* is said 
to be an invariant domain for the process Sh if and only if for any state Uh in A, the 
state Shiuh) is also in A. 

For scalar conservation equations the notions of invariant sets and invariant 
domains are closely related to the maximum principle, see Example 2.3 In the 


case of nonlinear systems, the notion of maximum principle does not apply and 
must be replaced by the notion of invariant domain. To the best of our knowl¬ 
edge, the definition of invariant sets for the Riemann problem was introduced in 
[Nishida( 1968)1, and the general theory of positive ly invariant regions was developed 
in [Chueh et al.(1977)Chueh, Conley, and Smoller] . Applications and extensions to 
numerical methods were developed in |Hoff(1979)[ HofF(1985)| and [Frid(200T)] . 

The invariant domain theory when m = 2 and d = 1 relies on the existence 
of global Riemann invariants; the best known examples are the hyperbolic systems 
of isentropic gas dynamics in Eulerian and Lagrangian form, see Example |2.4| and 
Lions et al.(1996)Lions, Perthame, and Souganidis . For results on general hyper- 
bolic systems, we refer to [Frid(2001)| , where a characterization of invariant domains 
for the Lax-Friedrichs scheme and some flux splitting schemes is given. In particular 
the existence of invariant domains is established for the above mentioned schemes for 
the compressible Euler equations in the general case m = d+2 (positive density, inter¬ 
nal energy, and minimum principle on the specific entropy), see [Frid(2001) Thm. 7 
and Thm. 8]. Similar results have been established for various finite volume schemes 
in two-space dimension for the Euler equations in [Perthame and Shu(1996)| Thm. 3]. 

The objective of this paper is to propose an explicit numerical method based on 


continuous finite elements to approximate (1.1) such that any convex invariant set of 


(1.1) is an invariant domain for the process generated by the said numerical method. 

To facilitate the reading of the paper we now illustrate the abstract notions of 
invariant sets and invariant domains with some examples. 

2.3. Example 1: scalar equations. Assume that m = 1 and d is arbitrary, 
i.e., (1.1) is a scalar conservation equation. Provided / G Lip(IR; K'^), any bounded 
interval is an admissible set for (1.1). Eor any Riemann data ul,uh, the maxi¬ 
mum speed of propagation in ( |^ is bounded by Ainax(wL, Uii;) := 
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where Umin = niin(ttL, ttmax = niax(uL,ttjf). If f is convex and is of class 
C'\ we have Ai„ax(uL,Ufl) = max(|n-/'(Mi)|, \n-f'{uii)\) if ri-f'{uL) < n-f'^un) and 
>ym&^{uL,UR) = n-{f {ul) - f {ur)) / {ul - Ur) Otherwise. Any interval [a, 6] C K is ad¬ 
missible and is an invariant set for ([TAjA.e., if Ur, ul € [a, 6 ], then a < u(n,UL, ur) < 
b for all times; this is the maximum principle. For any a < 6 S K, the interval [a, b] 
is an invariant domain for any maximum principle satisfying numerical scheme. Note 
that the maximum principle can be established for a large number of numerical meth¬ 
ods (whether monotone or not), see for example Crandall and Majda(198^. 


2.4. Example 2: p-system. The one-dimensional motion of an isentropic gas 
is modeled by the so-called p-system, and in Lagrangian coordinates the system is 
written as follows: 


( 2 . 10 ) 


dtv + dxu = 0 , 

dtu -\- dxp{v) = 0 , for (x,t) S MxK+. 


Here d = 1 and m = 2. The dependent variables are the velocity u and the specific 
volume V, i.e., the reciprocal of density. The mapping v i—>■ p{v) is the pressure and is 
assumed to be of class (K+; K) and to satisfy 

(2.11) p' < 0, 0 < p". 


A typical example is the so-called gamma-law, p{v) = rv~'^, where r 
Using the notation u = {v,u)^, any set A in (0,oo)xIR is admissible. 

Using the notation dp := \/—p'{s) ds, and assuming dp < 
has two families of global Riemann invariants: 

POO POO 

(2.12) wi( m) = tt-I- / dp, and W 2 {u) = u — / dp. 

J V J V 

Note that /]°°dp < oo if 7 > 1. If 7 = 1 we can use wi{u) = u — ^Jrlogv and 
W 2 {u) = u + y/r\ogv. Let a, & G M, then it can be shown that any set Aab G K+xK 
of the form 


> 0 and 7 > 1 . 
00 , the system 


(2.13) 


Aab ■= {u G K+xM I a < W 2 {u), Wi{u) < b} 


is an invariant set for the system ( 2 . 10 ) for 7 > 1 , see |Hoff(1985)| Exp. 3.5, p. 597] for 
a proof in the context of parabolic regularization, or use the results from [Young(2002)| 
for a direct proof. Moreover, Aab is an invariant domain for the Lax-Friedrichs scheme, 
see [Hoff(197^ Thm. 2.1] and |HofF(198^ Thm. 4.1]. 

Since in the rest of the paper the maximum wave speed is the only information 
we are going to need from the Riemann solution, we give the following result. 

Lemma 2.5. Let {vl,ul),{vr,ur) with vr,vl < 00. Then 


Ainax('^L;'I*i?) —^ 


^-p'{min{vL,VR)), if ur - ur > ^[vr - vr){p{vr) - p{vl)), 


\/—p'{v*), otherwise, 

where v* is the unique solution of (j){v) := /l(u) -I- fR(v) -\- ur — ur = 0 and 

'-\/(piv) - p{vz){vz - v), ifv<vz 

dp, ifv>vz- 


fziv) := 
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Upon setting := max{wi{uL),wi{uii)) and := mi-a{w 2 (uL),W 2 {uii)) we 

have also have < min{vL,vji,v*), i.e., XmaxiuLyUn) < \/—p'iv'^), where 


:= (7r)^ 


(7 - 


2 

(7-I) 


Proof. It is well know that the solution of the Riemann problem consists of three 
constant states ul, u* , and un connected by two waves: a 1-wave connects ul and 
u*, and a 2-wave connects u* and ur. Moreover, a vacuum forms if and only if 
lim„_j.+oo > 0, see |Young(2002)| for details. In the presence of vacuum the 
equation 4>(v) = 0 has no solutions and in this case we conventionally set v* := -l-oo 
and ^J—p'(v*) := 0. Note that since (j) is an increasing and concave up function with 
lim„_>o+ 4>iv) = — 00 , the solution v* is unique. We also have that the maximum 
speed of the exact solution is Aniax('itL,Mi^) = max(-\/—p '(wl), \/—p'(v*), \/ —p' (vr)). 
The only possibility for Amax(itL,MR) = \/—p'{v*) is if v* < min(wi,Uij), i.e., the 
solution contains two shock waves which is equivalent to (j){mhi(yL,vji)) > 0. Us¬ 
ing the definition of 4> we derive that \max{uL,UR) = \/—p'{v*) if and only if 
(l){mm{vL,VR)) = ul - Ur - a /(vl - vr){p{vr) - p{vl)) > 0. This finishes the proof 
of the first part of the lemma. 

The exact value of v* can be found using Newton’s method starting with a guess 
< V*. This guarantees that at each step of Newton’s method the estimated max¬ 
imum speed is an upper bound for the exact maximum speed. One can obtain such 
° by using the invariant domain property (2.13), 

' = 1(^2 thereby giving 


a guess V 
:= (v 


° := (z;°, vP) by w 


max 

1 


= wi{u'^) and 


i.e., we define the state 


= ( 7 ^)' 


(7 - 


(7-1) 


The invariant domain property guarantees that < v*. Hence, the result is estab¬ 
lished. □ 

Remark 2.1. Note that the estimate on Amax(ML, ur) given in Lemma 2.5 is valid 
whether vacuum is created or not in the Riemann solution. 

Remark 2.2. We only consider the case where both and ur., are not vacuum 
states in Lemma [2.5[ since the algorithm that we propose in this paper never produces 
vacuum states if vacuum is not present in the initial data. 

2.5. Example 3: Euler. Consider the compressible Euler equations 


m 


(2.14) 


dtC + Vific))=0, 


/(c) = ( m^^+pl 
'^{E+p) 


m I 

P 


where the independent variables are the density p, the momentum vector field m and 
the total energy E. The velocity vector field u is defined by u := mjp and the internal 
energy density e by e := E — The quantity p is the pressure. The symbol I 

denotes the identity matrix in Let s be the specific entropy of the system, and 
assume that —s{e,p~^) is strictly convex. It is known that 


(2.15) 


Ar '.= {(p, m, pE) I p > 0, e > 0, s > r} 
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is an invariant set for the Euler system for any r S K. It is shown in |Frid(200T)j Thm. 
7 and 8 ] that the set Aj. is convex and is an invariant domain for the Lax-Friedrichs 
scheme. 

Let n G 5'^“^(0,1) and let us formulate the Riemann problem (2.1) for the 


Euler equations. This problem was first described in the context of dimension split¬ 
ting schemes with d = 2 in [Chorin(1976) p. 526]. The general case is treated in 
[Colella(1990)[ p. 188], see also |Toro(2009)[ Chapter 4.8]. We make a change of ba¬ 
sis and introduce ti,..., so that {n, ti,..., td-i} forms an orthonormal basis 
of With this new basis we have m = where m := pu, u := M-n, 


■- piu-ti ,..., u-td-i) ■■= pu^. 

The projected equations are 

( P\ 


(2.16) dtc +d^{n-f{c)) =Q, 

c = 

m 

\E ) 

, n-f{c) = 

-I- p 
um^ 

\u{E+p)j 


Using the density p and the specihc entropy s as dependent variables for the pressure, 
p{p,s), the linearized Jacobian is 


/ u p 0^ ® \ 

p~^dpp u o''' p~^dsp 

0 0 ul 0 

\ 0 0 QT u / 


The eigenvalues are w, with multiplicity d, u + \jdpp{p, s), with multiplicity 1, and 
u — dpp{p, s), with multiplicity 1. One key observation is that the Jacobian does 
not depend on m-*-, see |Toro(20(i^ p. 150]. As a consequence the solution of the 
Riemann problem with data (c/,, C/j), is such that (p, u,p) is obtained as the solution 
to the one-dimensional Riemann problem 


(2.17) 



+p I = 0, 

u{£ +p)j 


m 


with e = £ — 

2p 


•'Z 11^2 


with data := {pL,mL-n,£L), := {pR^mR-n, £fi), where £z = Ez - 

Z G {L,R}. Moreover, for an ideal gas obeying the caloric equation of state p = 
(7 — l)pe, it can be shown (see [Toro(2009) p. 150]) that m,-^ is the solution of the 
transport problem dtm^ + dx{um) = 0. The bottom line of this argumentation is 
that the maximum wave speed in (|2.16|) is 


Amax(cL,Cfl) = max(lA;i (c2,c^)l, lA+(< 


■L^^R 


)!)• 


where A]"(c”,c)^) and Ag (c£,c)|) are the two extreme wave speeds in the Riemann 
problem (2.17) with data (c",c^). 

We now determine the values of A]"(c”,c)|) and Ag (c2,c^). We only consider 
the case where both states, cr and cr, are not vacuum states, since the algorithm 
that we are proposing in this paper never produces vacuum states if vacuum is not 
present in the initial data. That is, we assume PLtPr > 0 and PLtPr > 0. Then the 
local sound speed is given by az = \ where Z is either L or R. We introduce the 
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following notations Az := Bz '■= and the functions 

(2.18) (l){p) := f{p,L) + f{p,R)+ UR-UL 

ip-Pz)[^)^ Ap>pz, 


(2.19) 


f{p,Z) := 


where again Z is either L or R. It is shown in [Toro(2009)[ Chapter 4.3.1] that the 
function (j>{p) G C'^(M+;K) is monotone increasing and concave down. Observe that 
0(0) = Ur — ul — — 1^. Therefore, 0 has a unique positive root if and only if 

the non-vacuum condition 


( 2 . 20 ) 


2aL 2aR 
UR-UL< -7 + 


7 — 1 7 — 1 


holds, see |Toro(2009)t (4.40), p. 127]; we denote this root by p*, i.e., 0(p*) = 0 and 
p* can be found via Newton’s method. If (2.20) does not hold we set p* = 0. Then it 


can be shown that, whether there is formation of vacuum or not, we have 

1 

7-1-1 f p* -pl' 


( 2 . 21 ) 

( 2 . 22 ) 


-^1 — Ur — Ul 

+ Ofl [ 1 -f 


‘^1 \ PL / 

7 + 1 f P* -P r" 
27 V PR 


where z+ := max(0, z). 

Remark 2.3. (Fast algorithm) Note that if both 0 (pl) > 0 and 4>{pr) > 0, there 
is no need to compute p*, since in this case X)) {ur,ur) = ur — or and X^{ur, ur) = 
ur+clr, i.e., two rarefaction waves are present in the solution with a possible formation 
of vacuum. This observation is important since traditional techniques to compute p* 
may require a large number of iterations in this situation, see [Toro (200^ p. 128]. 


Note finally that there is no need to compute p* exactly since one needs only an upper 
bound on A^ax. A very fast algorithm, with guaranteed upper bound on A^ax up to 
any prescribed accuracy e of the type Amax + Amax + (1 + e)Amax, is described in 
Guermond and Popov(2015b)|. 


3. First order method. We describe in this section an explicit first-order finite 
element technique that, up to a CFL restriction, preserves all convex invariant sets of 


(1.1) that contain reasonable approximations of Uq. Although most of the arguments 


invoked in this section are quite standard and mimic Lax’s one-dimensional finite 
volume scheme, we are not aware of the existence of such a finite-element-based scheme 
in the literature. 


3.1. The finite element space. We want to approximate the solution of (1.1) 


with continuous finite elements. Let (Th)h>o be a shape-regular sequence of affine 
matching meshes. The elements in the mesh sequence are assumed to be generated 
from a finite number of reference elements denoted Ki ,..., For example, the 
mesh Th could be composed of a combination of triangles and parallelograms in two 
space dimensions (ci7 = 2 in this case); it could also be composed of a combination of 
tetrahedra, parallelepipeds, and triangular prisms in three space dimensions {vj = 3 
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in this case). The affine diffeomorphism mapping Kr to an arbitrary element K € Th 
is denoted Tk : Kr —K and its Jacobian matrix is denoted Jk, 1 < r < tu. We 
now introduce a set of reference Lagrange finite elements 5ir)}i<r<ro (the 

index r € {1:tu} will be omitted in the rest of the paper to alleviate the notation). 
Then we define the scalar-valued and vector-valued Lagrange finite element spaces 

(3.1) P{Th) = {vG C°(P;M) I v\koTk G P, VK g %}, P{Th) = [PiTh)r. 

where P is the reference polynomial space defined on K (note that the index r has 
been omitted). Denoting righ := dimP and denoting by {ai}ie{i:ra 3 h} the Lagrange 
nodes of K, we assume that the space P is such that 

(3.2) min v{a£) <v{x) < max v{ai), Vv G PjVx G K. 

l<^<nsh 

Denoting by Pi and Qi the set of multivariate polynomials of total and partial degree 
at most 1, respectively; the above assumption holds for P = Pi when K is a simplex 
and P = Qi when K is a parallelogram or a cuboid. This assumption holds also for 
first-order prismatic elements in three space dimensions. 

Let {ai}ie{i:/} be the collection of all the Lagrange nodes in the mesh 77i, and 
let be the corresponding global shape functions. Recall that 

forms a basis of P(Th) and ^i{aj) = Sij. The Lagrange interpolation operator in 
V{Th) is denoted IVh '■ C°{D) —s> T’(7h). Recall that n/i(u) = We 

denote by Si the support of ipi and by |5'i| the measure of Si, i G {!:/}. We also 
define S,j := Si n Sj the intersection of the two supports Si and Sj . Let P be a union 
of cells in Th', we define P{E) := {j G {!:/} | \Sj H E\ yf 0} the set that contains 
the indices of all the shape functions whose support on E is of nonzero measure. We 
are going to regularly invoke I{K) and P{Si) and the partition of unity property: 
J2iex{K) V5i(®) = 1 ab X G K. 

We define the operator C : P{Th) —^ so that C{vh) is the coordinate vector 
of Vh in the basis {^i}i(£{i:i}, i.e., Vh = Note that C(vh)i = Vh{ai). 

We are also going to use capital letters for the coordinate vectors to alleviate the 
notation; for instance we shall write V = C(vh) when the context is unambiguous. 
Note finally that the above assumptions on the mesh and the reference elements 
imply the following property: for a\\ x G K and all PT gTh, 

(3.3) min C{vh)t < Vh{x) < max C{vh)t, ^Vh G P{Th)- 

tei(if) tei(if) 

We define similarly the C : P{Th) —^ (M™)^, i.e., Vh = or equiva¬ 

lently C{vh)i = Vh{ai). 

Let AJ G be the consistent mass matrix with entries Jg ipi(x)(pj(x) dx, 

and let be the diagonal lumped mass matrix with entries 

(3.4) mi := / (pi(x)dx. 

■J Si 

The partition of unity property implies that mi = ) I ba:, i.e., the 

entries of are obtained by summing the rows of A4. 
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3.2. The scheme. Let Uho S ^{Ih) be a reasonable approximation of Mq (we 
shall be more precise in the following sections). Let n S N, r be the time step, 
be the current time, and let us set = t” + r. Let G 'Pi'Th) be the space 
approximation of u at time and set U" = Q{u^). We propose to compute by 

Ijn+l _ M" r 

(3.5) / V-(n^/(M))))(^,d^- V d,,UJ=0, 

where and the lumped mass matrix is used for the approximation 

of the time derivative. The coefficient dij is an artificial viscosity for the pair (i, j) 
that has yet to be clearly identihed. For the time being we assume that 

(3.6) dij ^0, if ^ ^ J; dij — dji^ and da .— ^ ) dji. 


Using that V-(n,i/(M”)) = J2j /(U” )-V(pj, the above equation simplihes into 
yn+l_ yn 


(3.7) 


- ^ /(up-Q, - Ujd,, = 0, 

jex{s,) 


where the coefficients G are defined by 


(3.8) 


Cij = / ipi\7(pj dx, 


ID 


Remark 3.1. (Conservation) The definition da := J2i^j^x{s ) ~^3i impbes that 
SjGi(Si) which in turn implies conservation, i.e., dx = dx + 

-{Ilhfiu]^)) dx. Note also that the symmetry assumption in (3.6) implies da := 
'^i^j^i{Si) which is often easier to compute. 

3.3. The convex combination argument. We motivate the choice of the 
artihcial viscosity coefficients dij in this section. Observing that the partition of 
unity property J2j^x{Si) ~ ^ (3.6) imply conservation, i.e., 

(3.9) Cy = 0, dij = 0. 

i6X(Si) iGl(Si) 


we re-write (|3.7|) as follows: 
(3.10) 


i|"+i _ II" 

= - E (/(u”)-/(ur))-Q,+d.,(u” + ur). 

1 GI(S.) 


Using again conservation, i.e., da = — J2i^j^x{S ) finally arrive at 


(3.11) 


U 


n+1 


ur(i- E 


2Tdi 


E ^xdij jj"+i 
m.. b 


ii^jexiSi) ii^jexis,) 

where we have introduced the auxiliary quantities 

:= + UD - (/(UJ) - /(UD)-^ 


(3.12) 
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A first key observation is that (3.11 1 is a convex combination provided r is small 


enough. A second key observation at this point is that upon setting := Cjj/||cij||£ 2 , 
t = ||( 

the following definition for the viscosity coefficients dtj 


is exactly of the form it(t, riij, Ui, Uj) as defined in ( 2 . 6 ) with a fake time 


\p/2dij. The CFL condition tAi„ax(^iy > mlj Mfl) < 5 in Lemma 


2.1 


motivates 


(3.13) 


dij .— max(Aniax(^ij ; , U j 


|^2,An,ax(n,*,U”,Ur)||c,y||,2), 


where recall that Amax(^y) U^, Uj) is defined in the assumption (2.5). 

Remark 3.2. (Symmetry) If either or aj is an interior node in the mesh, one 
integration by parts implies that Cij = —Cji, which in turn implies Xmaxinij, Ui, Uj) = 
Aniax(^ji 5 bij , U j) . In conclusion Amax(^2j5 — Aniax(^_;z) bij , U j ) II II ^2 

if either ai or aj is an interior node. 

Remark 3.3. (Upwinding) Note that in the scalar one-dimensional case when the 
flux / is linear, (|3.5|) gives the usual upwinding first-order method. 


4. Stability analysis. We analyze the stability properties of the scheme (3.5) 


with the viscosity defined in (3.13). 


4.1. Invariant domain property. Upon defining hx '■= diam(A'), the global 
maximum mesh size is denoted h = max/^-gTy Hk- The local minimum mesh size, 
for any K GTh is defined as follows: 


(4.1) 


hx ■— 


max,^jgx(if) ||V(/?*|l£=o(s,^.) ^ 


and the global minimum mesh size is h := uAnx^Th ttK- Due to the shape regularity 
assumption, the quantities hx ^^.d hx are uniformly equivalent, but it will turn 
out that using hx instead of hx gives a sharper estimate of the CFL number. Let 


rish := card(I(Ar)) and let us define dx '■= 


"sh-l 


. Note that 


(4.2) 


0 < i^min := min min dx < +00, 

(J~h)h>o 


since there are at most w reference elements defining the mesh sequence. We also 
introduce the mesh-dependent quantities 

(4.3) /Tinin := min min / (pi{x)dx, ^max := max max / (pi(x)dx. 
KeThiGliK) \K\ Jx K(iThr£l(K)\K\Jx 

Note that pmin = hmax = for meshes uniquely composed of simplices and 

Akinin = Mmax = for meshes uniquely composed of parallelograms and cuboids. We 
now prove the main result of the paper. 


Theorem 4.1. Let A C A be an invariant set for (1.1) in the sense of Defini- 
tion \2.fA Assume that A is convex and 

(4.4) Amax(A) := max max Ainax(w, t6x, m_r) < 00, 

neS‘^-i(0,l) ULjUneA 

where S'‘^“^( 0 , 1 ) is the unit sphere in Assume that Uho S A and r is such that 


(4.5) 


2r 


'^max (A) n 

max 

h hmin'dn 


< 1 . 


Then 
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(i) A is an invariant domain for the solution process i —> for all n > 0. 

(ii) Given n > 0 and i € {1:J}, let B C A be a convex invariant set such that 
Ur G B for all I G I{Si), then G B. 

Proof. We prove the statement 0 by induction. Assume that G A for some 
n > 0; we are going to prove that G A. Note that Uho G A by assumption. Let 
i G {1:/} and consider the update (3.5) rewritten in the form (3.11). Observe that 


772+1 


\cij\\p/2dij. The definition 


upon defining := Cy7||cij||£2, the quantity U+ defined in (3.12) is exactly of 
the form it(t, U”, U") as defined in (2.6) with the flux f riij and the fake time 

t = 


(4.6) 


dij P Aj] 




I 


is the CFL condition for the conclusions of Lemma 12.11 to hold with fake time t = 

— —n+l 

\\cii 11/2 /2c?,;i. Since A is a convex invariant set we have Uj := u{t, riij, U", U") G A 


for all j G I (Si). Let us now prove that (3.11) is indeed a convex combination by 
proving that 1 - = 1 + > 0. Note first that 

— I II II £2 (pi dx ^ h / ip^dxtfh /Tniax I I. 

JSii JSij 


The definition of da implies that 


7 ^ Amax(A) 

-dti < - - - n 


max E i^bi< 


'^max(-^) n 

h dr, 




The using that /.tmin|<5'i| < mi, we infer that 


_2t— < 2r 
mi 


Amax(A) 11 

n 


M: 


ndr, 


< 1 , 


mm ^ min 


in+l 


which proves the result owing to the CFL assumption (4.5). Hence (3. Ill) defines UJ 

as a convex combination between U" and the collection of states {Uj^ }jei{Si)- The 

convexity of A implies that G A, since U" G A by assumption and we have 

-? 2 “|“ 1 

established above that Uj^ G A for all j G X{Si). The space approximation being 
piecewise linear, a convexity argument implies again that G A, which proves the 
induction assumption. 

Note in passing that we have also proved the following local invariance property: 
given any convex invariant set 5 C A that contains {U"}/gx(Si)i in B, 

i.e., the local statement (|^ holds. This completes the proof. □ 


Remark f.l. The arguments invoking the convex combination (3.11) and the one¬ 
dimensional Riemann averages ( 3.12[) are s imilar in spirit to those used in the proof 
of Theorem 3 in Perthame and Shu(1996)| . 

4.2. Discrete entropy inequality. We now derive a local entropy inequality. 


Theorem 4.2. Let A G A he a convex invariant set for (1.1). Let {r],q) be an 
entropy pair for (1.1). Assume that (2.7) holds for any Riemann data {uL,Ufj) in 
A, and any n G 1). Assume also that (4.4) and (4.5) hold, then we have the 

following for any n > 0 and any i G {1:7}-' 


(4.7) !+(r;(U+')-77(Ur))+ / ViIVnq{ul))ip,dx 

T Jd 


E ^b^(U”) < 0. 
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Proof. Let (?7, <7) be an entropy pair for the system (1.1). Let i S {!:/}, then 

;hat 

‘Irdij /yrn+l, 

^(U„- ). 


recalling (3.11), the CFL condition and the convexity of rj imply that 

2 rd 


i?(ur')<(i- ^ 

Owing to Lemma [2.2| we have 


^)t7(ur)+ ^ 




?7(ur/^) < 5(r?(ur) +?7(up) -q(un-nzi). 

with t = \\cij\\i 2 /2dij] hence, 


- r;(Ur)) < 5: 2d,,(iy(U:;+^) - r?(U”)) 

i¥=3eI{Si) 

< Y. ^d(^(U ■) - vm)) - ||o,||,.(q(UJ)-n,, - g(U").n.,). 


The conclusion follows from the definitions of riij, c^j and dy. 


Remark 4-2. One recovers the equation (3.5) from \A.1\ with ri[v) = v. Note also 
that (4.7) gives the global entropy inequality 771^77(1)"“''^) < X]i<i</ 


Remark 4 . 3 . The meaning of the entropy inequality (2.7) might be somewhat 
ambiguous in some cases, especially when m is a measure. Since it is only the in¬ 
equality (2.8) that is really needed in the proof of Theorem 4.2 we could replace 
the assumption (|2.7|) by (2.8). This would avoid having to invoke measure solutions 


since u{t, n, uh) should always be finite for the Riemann problem (2.1) to have a 
reasonable (physical) meaning. 


4.2.1. Cell-based vs. edge-based viscosity. In the formulation (3.5) the 
term X]jGi(Si)models some edge-based dissipation, i.e., dij is a dissipation 
coefficient associated with the pair of degrees of freedom of indices {i,j). This for¬ 
mulation is related in spirit to that of local extremum diminishing (LED) schemes 
developed for scalar conservation equations in [Kuzmin and Turek(2002)1 Eq. (32)- 
(33)], see also |Jameson(1995)[ §2.1]. It is however a bit difficult to understand that 
we are modeling some artificial dissipation by just staring at ( |3.5| ). 

We now propose an alternative point of view using a cell-based viscosity. The 
traditional way to introduce dissipation in the finite element world consists of invoking 
the weak form of the Laplacian operator —V-(i7V'0). Eor instance, assuming that the 
viscosity field v is piecewise constant over each mesh cell K € Th, we write: 


(4.8) 


/ -V-{vV'ilj)ipidx = 
d-D KcSi ^ 


vk I Vif-Vipidx. 

' K 


Unfortunately, it has been shown in Guermond and Nazarov(2013)| that the bilinear 
form VV'-Vi^da: is not robust with respect to the shape of the cells. 

More specifically, the convex combination argument, which is essential to prove the 
maximum principle for scalar conservation equations in arbitrary space dimension 
with continuous finite elements, can be made to work only if Vipi-Vipj da; < 0 for 
all pairs of shape functions, Lpi, (fj, with common support of nonzero measure. This 
is the well-known acute angle condition assumption, which a priori excludes a lot of 
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meshes in particular in three space dimensions. To avoid this difficulty, it is proposed 
in Guermond and Nazarov(2013)| to replace (4.8) by J2kcS where 


(4.9) 


bK{(pj,(pi) = {\K\ iti=j, i,jGl{K), 

0 or j ^I{K). 


The essential properties of bx can be summarized as follows: 

Lemma 4.3. There is c> Q depending only on the collection {{K^, P^, L]r)}i<r<CT 
and the shape-regularity, such that the following identities hold for all K G Ph and all 

Uh,Vh€ P{Th): 

(4.10) bK{(pt,‘Pj) = bK{(pj,(pi), bxiTi, X! Tj)=0. 

(4.11) bK{uxVk) = dK\K\ ^ ^ (U,-U,)(V,-V,) 

(4.12) bK{uh,Uh) > ch^||Vuft,||^2(^). 

For instance, when iV is a simplex and K is the regular simplex, i.e., all the edges 
are of unit length, it can be shown that bxiTiTTi) = ^ Jk da; and 

bK{(pj,<Pi) = /k J]f(Vv?j)-Jl-(V(^i)da; for j ^ i, with k = Kl + g). Note also 

that bxi^PjT^Pi) ~ hf^ fj^(V(pj)-(V(pi)dx if iV is a regular simplex, thereby showing 
the connection between bx and the more familiar bilinear form associated with the 
Laplacian. One key argument from Guermond and Nazarov(2013)| is the recognition 


that the bilinear form defined in (4.9) has all the good characteristics of the Laplacian- 


based diffusion (see Lemma 4.3) and makes the convex combination argument to work 
independently of the space dimension and the shape-regularity of the mesh family. 


Hence, instead of (|3.5[), we could also compute by 


(4.13) 


U 


n+1 


u( 


+ [ V-{U!,fiK))^,dx+ Y. Y. ^7>^KiTJ,T^)=0 


KeTh jGi(if) 


where {^')(-}iCGTh is a piecewise constant artificial viscosity scalar field. 
Theorem 4.4. Let Wx}K&n be defined by 


(4.14) 


i^x ■= 


max 


-^max {riij, Ui, Uj)|jcy 11^2 


Then the conclusions of Theorem \41\ and Theorem 14-^ hold under the assumptions 
(4.4) and (4.5) and with the solution process itj) — n>0, defined by (4.13). 

Proof Let us denote dij := —J2 xgS then (4.13) can be recast as 

follows: 


U 


n+1 


ui 


- E 

j&ps,} 


U”dy = 0, 


which in turn implies that 

J”( 


U 


n+1 


= u” 1- 


E 




2Tdi 


rui 




2 t dij 
mi 
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where we have introduced the auxiliary quantities 


:= hu] + UD - (n,,-/(Up - 

^ ^cli 


'IJ 




yn+l 


Here again is of the form u{t, Ui, Uj) as defined in (2.6) with the fake time 


t = \\cij\\g 2 / 2 dij, hence we need to make sure that Ai„ax(’^y; ||^2/2dy < ^ 

to preserve the invariant domain property. Recalling that dij has been defined by 
dij := Xma-Kinij,UL,uii)\\cij\\i 2 (see ( |3.6[ )), the above condition reduces to showing 
that dij < dij. The definitions of vk and dij implies that 


dij = - ^KbK{^j,^z)> - 






KGSi. 


> — ^ _ _ 

Etcs,, 


iCGS,, ^rcS.,- ^T{‘Pj,V’i) 

--bK{ifj,tpi) = dij, 




whence the desired result. We now prove that 1 —under the CFL 


condition (4.5). From the proof of Theorem 
hence 


4.1 


we have dij < An 


C {A)h- 

/^max 


U h 


I'K < 


'^max(^)/^ n 


max 


\Su\ 


k^lGl{K) X^TC^fci ^T{^kf 


which in turn implies that 


dij ^ 


•^niax {A)ii 

n 


h 


E 

KcS,j 


-bK{(fi,(pj) 


max — - 
ki^lGl(K) EtcSh 


\Su\ 

-bT{vk,Vi) 


Recalling the definition of bT{{^Pk,^i) we have X^TcSti ‘i’*) — '*^min|T| = 

^min \Ski\; hence 


'^max(-^)/^ n 


' •max \/^^iiiax \ ^ 7 / \ 

ij< - -^—T - 2^ -bK{Vi,Vi) = 

k' min n 


K(ZSi. 


'^max(-^)M n 

79 • h 

mm ' 


E ^^1^1- 


KcSi. 


Finally we have 

da .— 2 ^ dij ^ 


^maxiA') fl 

n 

79 . /E 




i#jGl(Si) KcSij 


i") ■ h 

^minLF 


This means that the bound on —da is the same as that on —da in the proof of 
Theorem 4.1 This concludes the proof. □ 


5. Numerical illustrations. We illustrate in this section the method described 
in the paper, i.e., (3.5)-(3.13), and discuss possible variants. 


5.1. Invariant domain property and convergence issues. We give in this 
section a counter-example showing that a method that is formally first-order consistent 
and satisfies the invariant domain property may not necessarily be convergent. 

To illustrate or point, let us focus our attention on scalar conservation equations 
and let us consider an algebraic approach that is sometimes used in the literature, see 
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e.g., [Kuzmin et al.(2005)Kuzmin, Ldhner, and Turek p. 163], [Kuzmin and Turek(2002)1 | 
Eq. (32)-(33)]. Instead of constructing a convex combination involving (entropy sat¬ 
isfying) intermediate states like in (3.111, we re-write (3.10) as follows: 


(5.1) 


U«+i_u- 


= - E E 


i¥^jex{Si) 


jex{s,) 


Or, equivalently 

U"+^ - U 

(5.2) 


- = - E 

^^j&X{Si) 


fm-fm 


U^-U^ 


Q ,( u -- ur )+ ^ 


jexis.) 


Let us set kij := ^ -Cij, (with kij := 0 if U” = U”), then 


(5.3) ur^ = ur(l-^ E i-k^,+d^,))+ E 

\ II L'i / II L't 


i^jex(Si) 


i¥^jex{Si) 


Let us finally set 

(5.4) dij := max{0, kij, kji), i^j, and du := — d, 


t] ■ 


This choice implies that —kij 


U 


n+1 


dij > 0 for all i 


G {1:1V}, j G 2{Si). As a result, 


G conv{U”, j G T(S'i)} under the appropriate CFL condition; hence, the solu¬ 
tion process i— described above in (5.3)-(5.4) satisfies the maximum princi¬ 
ple. Although, this technique looks reasonable a priori, it turns out that it is not dif¬ 


fusive enough to handle general fluxes as discussed in Guermond and Popov(2015a) 
§3.3]. The convergence result established in Guermond and Popov(2015a)| requires 


an estimation of the wave speed that is more accurate than just the average speed 
" ^ ’ which is invoked in the above definition. This definition of the wave 

speed is correct in shocks, i.e., if the Riemann problem with data (U^, U_,) is a sim¬ 
ple shock; but it may not be sufficient if the Riemann solution is an expansion or a 
composite wave, which is likely to be the case if f is not convex. 

We now illustrate numerically the observation made above. We consider the so- 
called KPP problem proposed in Kurganov et al. (2007) Kurganov, Petrova, and Popov] .| 
It is a two-dimensional scalar conservation equation with a non-convex flux: 


(5.5) 


dtu + V • f{u) = 0, u(x, 0) = uo{x) = 


147r 

4 


4 ’ 


if y/x'^ + < I, 

otherwise. 


where f{u) = (sin u, cos u). This is a challenging test case for many high-order numeri¬ 
cal schemes because the solution has a two-dimensional composite wave structure. For 


example, it has been shown in Kurganov et al.(2007)Kurganov, Petrova, and Popov 


that some central-upwind schemes based on WEN05, Minmod 2 and SuperBee re¬ 
constructions converge to non-entropic solutions. 

The computational domain [—2,2] x [—2.5,1.5] is triangulated using non-uniform 
meshes and the solution is approximated up to t = 1 using continuous Pi finite ele¬ 
ments (29871 nodes, 59100 triangles). The time stepping is done with SSP RK3. The 


solution shown in the left panel of Figure 5.1 is obtained using (5.4) for the definition 
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of dij. The numerical solution produces very sharp, non-oscillating, entropy violating 
shocks, the reason being that the artificial viscosity is not large enough. Note that the 
solution is maximum principle satisfying (the local maximum principle is satisfied at 
every grid point and every time step) and no spurious oscillations are visible. The nu¬ 
merical process converges to a nice-looking (wrong) piecewise smooth weak solution. 
The numerical solution shown in the right panel of Figure [5T] is obtained by using our 
definition of , ( 3.13) (note in passing that the results obtained with (4.13l-(4.14| 
together with (3.13) are indistinguishable from this solution). The expected helicoidal 
composite wave is clearly visible; this is the unique entropy satisfying solution. 



Fig. 5.1. KPP solution with continuous Pi elements (29871 nodes, 59100 triangles). Left: 
entropy violating solution using ( |3.5[ |-( |5^ ,' Right: entropy satisfying solution using 


In conclusion, the above counter-example shows that satisfying the invariant do¬ 
main property/maximum principle does not imply convergence, even for a first-order 
method. It is also essential that the method satisfies local entropy inequalities to be 


convergent; this is the case of our method (3.5)-(3.13) (see Theorem 4.2), but it is not 


the case of the algebraic method (5.3)-(5.4). 


Remark 5.1. The reader should be aware that we are citing Kuzmin et al.(2005)Kuzmin, Lohner, and Turek ,| 
p. 163], [Kuzmin and Turek(2002)[ Eq. (32)-(33)] a little bit out of context. The 
scheme as originally presented in the above references was only meant to solve the 
linear transport equation, and as such it is a perfectly good method. Problems arise 


with (5.4) only when one extends the methodology to nonlinear nonconvex fluxes, as 


we did in (5.2). 


5.2. Special meshes. The construction of the intermediate states in ( 3.10| ) is 
not unique. For instance we can extend a construction used by |Hoff(1979) Cor. 
1] in one space dimension for the p-system. Let us assume that i S {1,...,1V} is 
such that every j G X{Si) \ {i}, there is a unique <Ji{j) G I{Si) \ {i,j} such that 
Cij := Jg (piV(j)jdx = — dx =: This property holds in one 

space dimension for any mesh if is an interior node. It holds in higher space 
dimension provided the mesh has symmetry properties and is an interior node; for 
instance it holds if the mesh is centrosymmetric, i.e., the support of (j)i is symmetric 
with respect to the node ai for any i G {1,...,1V}. Then we can re-write (3.7) as 
follows: 


(5.6) rrii 


U 


n+1 


ui 


= d,.ur- 5] (/(u")-/(u" (,.)))• 

leJiSi) 


i d-d. 
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where the set J{Si) C T{Si) is such that ai : J{Si) — (Ji{J{Si)) is bijective and 
J{Si) U (Ji{J{Si)) = I{Si) \ {i}. Then upon recalling that da := + 

d'iaiU)), 


(5.7) Ur =Ur 1- 


E 

m. 




■ ^iCTiU)) 


E 


T{di. 


Tj (j)) 77”+l 

'■’ij > 


_71 ~ 1 ” 1 

where we have defined the intermediate state Uj. by 


(5.8) UT = 


The state is of the form u{t, U") := /“” u{nij, U”)(a:, t) dx 


Hcri{j) 


-U 


• O') 


UY 


dij di(j^(^j'j 


u;‘-(/(U 7 )-/(u-(^.))) 


^U) 


where aL = — 


. , . , OiR= , ,y - and t := , , provided 


iU) 


(5.9) > (Ar)-(n,,, U;)||c„-||,., Vj S J{S.), 

(5.10) d,, > (A+)+(n,„ Up|lc,,||,., Vj e J(^,), 

where we defined x'^ = max(a;, 0) and x~ = — min(a;,0). A sufficient condition that 
implies both the above inequalities and is independent of the choice of the set Yfi{Si) 
is 


(5.11) 


min(dy,dio..(j)) > A 

n 




u( 


”ij\\P ! 


j G J{S,). 


Note that the above argument holds only if ai is an interior node satisfying the 
symmetry property If this is not the case, then we can always use the 

lower bound ( [T6| ), i.e., d^j > Ainax(w*i, Uf, U")||c,; 

In conclusion the diffusion matrix 


II II 

"y! ‘-’i j '-'j JU'-tj I 
{dij)i<ij<N 


can be constructed as follows: 


(1) For every node i satisfying the symmetry property 
G f7(5*j), we define dij — dio-^^j) — Aniax 


= -CiaiU) 

3 G jpij, we oenne cty = a^ai(j) = ^max(,«ij, uyj||cij ||£ 2 ; (2) For every 

(rther index i not satisfying the symmetry property mentioned above, we define 
dij = Xmaxinij, U", U")||cy 11 ^ 2 ; (3) We construct the diffusion matrix by setting dij := 
max(dy, dji) for j ^ i and da := — '^i^j^x{Si) This construction guarantees con¬ 
servation, i.e., J2i£X(s-) ~ ^ first-order consistency, i.e., J2j£X(Si) = 0. 

Remark 5.2. Quite surprisingly, in the case of scalar linear transport the above 
construction and the construction done in (3.3 (see definition (3.13)) give the same 
scheme (i.e., the same CFL). 


5.3. Invariant domain property vs. monotonicity. We show in this sec¬ 
tion that the invariance property and what is usually understood in the literature 
as monotonicity are two different concepts and just looking at monotonicity may be 
misleading. 


5.4. p-system. We consider the p-system and solve the Riemann problem cor¬ 
responding to the initial data {vl,ul) = (1,0), = (2^^, ;:^^). The compu¬ 

tational domain is the segment [0,1] and the separation between the left and right 
states is set at Xq = 0.75. The solution is a single rarefaction wave from the first 
family (i.e., Wi(vl,ul) = Wi{vn,UR)): 


[ 1 if 3^ < -1 

(5.12) v{x,t) = < 

I otherwise 
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(5.13) 



if < -1 

if _i < < _2- 

otherwise 


. 7 + 1 


This case is such that {v*,u*) = hence the second wave corresponding to 

the eigenvalues is not present. We use continuous piecewise linear finite elements 
with the algorithm ( |3.5| )-( 3T^ . The time stepping is done with the SSP RK3 tech¬ 
nique. We show the profile of w at t = 0.75 in Figure 5.2 for meshes composed of 
10^, 2x 10^, 4x 10^, 10^, 2x 10"*, 4x 10^, 10®, 2x10® cells. We observe that the profile 




Fig. 5.2. Left: v-profile for the p-system at t = 0.75, 10® grid points. Right: close up view of 
the V-profile for various grid sizes: 10®, 2x 10®, 4x 10®, 10"*, 2x 10^, 4x 10^, 10® grid points. 


is not monotone. There is an overshoot at the right of the foot of the (left-going) 
wave. Actually this overshoot does not violate the invariant domain property; we 
have verified numerically that, at every time step and for every grid point in each 
mesh, the numerical solution is in the smallest invariant domain of type (2.13) that 
contains the piecewise linear approximation of the initial data. This result seems a 
bit surprising, but it is perfectly compatible with Theorem |4.1| Since the numerical 
solution cannot stay on the exact rarefaction wave (green line connecting Ul and Ul 
in Figure 5.3), the second wave reappears in the form of an overshoot at the end of 
the rarefaction wave (see right panel of the Figure 5.2). 


U Ur 



V 


Fig. 5.3. The overshooting mechanism for a single rarefaction wave in the phase space for the 
p-system. Initial data in black; additional points after one time step in red; after two time steps in 
blue. Observe the position o/Ug. 

Let (Ul ■.., Ul, Ulj ..., U_r) be the initial sequence of degrees of freedom. Af¬ 
ter one time step two additional points appear in the phase space, denoted on Fig- 
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ure 


5.3 


by Uj and 1 ) 2 - Because of the invariant domain property, these points 
are under the rarefaction wave. Then the sequence of degrees of freedom at time 
t = T is (Ul ..., Ul, U}, U 2 , Ufj..., Ufl). Six additional points U^,. ■., Ug appear 
after two time steps and the sequence of degrees of freedom at time t = 2t is 
(Ul ..., Ul, U^, ..., Ug, U_R ..., Ui^). The point Ug is the one whose ^-component 
may overshoot because the exact solution of the Riemann problem with the left state 
U 2 and the right state U/j is composed of two rarefaction waves and the maximum 
value of V on these rarefactions is necessarily larger than vr (see red line in Fig¬ 
ure 5.31. Note that this is not a Gibbs phenomenon at all; in particular the amplitude 


of the overshoot decreases as the mesh is refined as shown in the close up view in 
the right panel of the Figure |5.2[ This phenomenon is actually very common in 
numerical simulations of hyperbolic systems but is rarely discussed; it is sometimes 
called ’’start up error” in the literature, see for example the comments on page 592 
Kurganov and Tadmor(2002)| and the comments at the bottom of page 1005 in 

The (relative) L^-norm of the error on both v and u at 
The method converges with an order close to 0.9. 


Liska and Wendroff(2003) 


t = 0.75 is shown in Table 


5.1 


\/h 

V 

rate 

u 

rate 

HP 

1.8632(-2) 

- 

7.2261(-3) 


2 x10^ 

1.0350(-2) 

0.85 

3.9239^3) 

0.88 

4x103 

5.6769(-3) 

0.87 

2.1173(-3) 

0.89 

OF 

2.5318(-3) 

0.88 

9.2888^4) 

0.90 

2 x10^ 

1.3644(-3) 

0.89 

4.9541(-4) 

0.91 

~AxW~ 

7.3151(-4) 

0.90 

2.6319^4) 

0.91 

1x103 

2.9695(-4) 

0.98 

1.1352(-4) 

0.92 

2x103 

1.5838(-4) 

0.91 

5.9869^5) 

0.92 


Table 5.1 

Convergence rates for the p-system 


5 . 5 . Euler in ID (Leblanc shocktube). We consider now the compressible 
Euler equations. We solve the Riemann problem also known in the literature as the 
Leblanc Shocktube. The data are as follows: 7 = § and 

Pl = 1 . 000 , ul = 0 . 0 , pl = 0.1 
Pn = 0 . 001 , Ur = 0 . 0 , pr = 10 “^®. 

The structure of the solution is standard; it consists of a rarefaction wave moving to 
the left, a contact discontinuity in the middle and a shock moving to the right. The 


density prohle is monotone. We solve this problem with the algorithm (3.5|-(3.13| 
using piecewise linear finite elements. The density profile computed with 50,000, 
100,000, 200,000, 400,000 and 800,000 grid points is shown in the left panel of Fig¬ 
ure |5.4[ The right panel in the figure shows a close up view of the region at the foot 
of the expansion wave. Of course the scheme does not have any problem with the 
positivity of the density and the internal energy, but we observe that the numerical 
profile is not monotone; there is a small dip at the foot of the expansion. There is 
nothing wrong here, since, for each mesh, the numerical solution is guaranteed by 
Theorem 14.11 to be in the smallest convex invariant set that contains the Riemann 
data. This phenomenon is similar to what has been observed for the p-system in the 
previous section. This example shows again that the invariant domain property is a 
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Fig. 5.4. Left: Density profile for the Leblanc Shocktube at t = 0.1. Right: close up view of the 
density profile at the foot of the rarefaction wave. 


different concept than monotonicity, and just looking at monotonicity is not enough 
to understand hyperbolic systems. 


6 . Concluding remarks. We have proposed a numerical method to solve hy¬ 
perbolic systems using continuous finite elements and forward Euler time stepping. 
The properties of the method are based on the introduction of an artificial dissipa¬ 
tion that is defined so that any convex invariant sets is an invariant domain for the 
method. The main result of the paper are Theorem 4.1 and Theorem |4. 2 1 The method 


is formally first-order accurate with respect to space and can be made higher-order 
with respect to the time step by using any explicit Strong Stability Preserving time 
stepping technique. Although, the argumentation of the proof of Theorem |4.1| re¬ 
lies on the notion of Riemann problems, the algorithm does not require to solve any 
Riemann problem. The only information needed is an upper bound on the local max¬ 
imum speed. Our next objective is to work on a generalization of the FCT technique 
(see Kuzmin et al.(2005)Kuzmin, Lohner, and Turek| ) to make the method at least 
formally second-order accurate in space and still be domain invariant. 
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